########################################################################
#   Project: The Press-Safety Paradox of Democracies in Media Systems, #
#            Regime-type Durability and Journalist Killings            #
#   Author: Jonathan A. Solis                                          #
#   Task: Replication Materials for Foreign Policy Analysis article -  #
#         Survival Analysis (manuscript and appendix)                  #
#   R version: 3.6.0                                                   #
#   Last updated: 5/31/2019                                            #
########################################################################

#----------------------------------------------------------------------------------------------#
# Replication Instructions: Use the Stata do.file Solis_fpa_may2019.do to create the following #
# .csv documents:                                                                              #
#                                                                                              #
# `survival_one.csv'                                                                           #
# `survival_two.csv'                                                                           #
# `survival_three.csv'                                                                         #
#                                                                                              #
# NOTE: I initially ran the survival analysis in R 3.5.1. Now that I've updated to R 3.6.0     #
# some of the results differ after the hundredths place in some analyses. However, the results #
# and inferences remain the same.                                                              #                                                               #
#----------------------------------------------------------------------------------------------#

#####################
# Survival Analysis #
#####################

#Create data for Table , Manuscript

rm(list=ls())

setwd("working director")

library(foreign)
library(survival)
library(norm)
library(ggplot2)
library(survminer)
library(GGally)
library(gmodels)

#                                                                                             #
# Data for Table 4, Manuscript:                                                               #
# Includes: Cox Proportional Hazards models 11-14, Therneau-Grambsch (TG) Tests, and          #
# calculations of Hazard rates.                                                               #
#                                                                                             #

#Load data
sa.one.full<-read.csv("survival_one.csv") 
names(sa.one.full)

#Stratify sample

#Autocracy
sa.one.auto <- subset(sa.one.full, durable2 == 0)
#Anocracy
sa.one.ano <- subset(sa.one.full, durable2 == 1)
#Democracy
sa.one.demo <- subset(sa.one.full, durable2 == 2)

############################
#####Full Sample############
############################

#Estimate model
cox.full <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.full, 
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.full

#                    coef  se(coef)       se2     Chisq   DF       p
#seq_ln          -0.23807   0.05927   0.05452  16.13429  1.0 5.9e-05
#polity2          0.01605   0.02008   0.01848   0.63891  1.0  0.4241
#public_cor       0.63357   0.40088   0.33114   2.49777  1.0  0.1140
#physical_vd     -3.50215   0.56520   0.48551  38.39400  1.0 5.8e-10
#express_vd       2.45502   0.62242   0.54669  15.55778  1.0 8.0e-05
#intensity2       0.83203   0.09106   0.08263  83.49428  1.0 < 2e-16
#info             0.01629   0.00574   0.00438   8.05830  1.0  0.0045
#pop_ln           0.43024   0.06359   0.04571  45.77353  1.0 1.3e-11
#frailty(ccode)                               167.31368 70.5 7.6e-10

#Iterations: 7 outer, 36 Newton-Raphson
#Variance of random effect= 0.637   I-likelihood = -1528.1 
#Degrees of freedom for terms=  0.8  0.8  0.7  0.7  0.8  0.8  0.6  0.5 70.5 
#Likelihood ratio test=947  on 76.3 df, p=0
#n=3586 (706 observations deleted due to missingness)

#T-G Test for proportionality 
cox.zph(cox.full, transform='log')

#                 rho  chisq      p
#seq_ln       0.01993  0.256 0.6127
#polity2     -0.04316  1.206 0.2722
#public_cor   0.02760  0.708 0.4001
#physical_vd -0.06890  4.015 0.0451
#express_vd   0.05005  1.926 0.1652
#intensity2  -0.07734  4.288 0.0384
#info         0.00883  0.091 0.7629
#pop_ln      -0.01882  0.412 0.5207
#GLOBAL            NA 13.298 0.1020

#Graph Scoenfeld Residuals
ggcoxdiagnostics(cox.full, type = "schoenfeld",
                 linear.predictions = FALSE,ggtheme = theme_bw())

#Result: physical_vd and intensity2 violate proportionality assumption.
#Diagnostics: Overall, model should be good because Global (p >.05). 

#ProportionalHazard
coef <- summary(cox.full)$coefficients[1, 1]
coef
exp(coef)
#0.7907903
1-0.7907903
#0.2092097; ~21%

##########################
#####Autocracy############
##########################

#Estimate model
cox.auto <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.auto,      
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.auto

#                  coef se(coef)     se2   Chisq  DF       p
#seq_ln         -0.5087   0.2018  0.1718  6.3571 1.0   0.012
#polity2         0.1422   0.2323  0.2035  0.3746 1.0   0.541
#public_cor     -0.4499   1.2867  1.0472  0.1222 1.0   0.727
##physical_vd   -2.2941   1.4752  1.2978  2.4182 1.0   0.120
#express_vd      0.3858   1.6088  1.4134  0.0575 1.0   0.810
#intensity2      1.3394   0.2603  0.2368 26.4782 1.0 2.7e-07
#info            0.0296   0.0167  0.0136  3.1320 1.0   0.077
#pop_ln          0.2848   0.1524  0.1273  3.4904 1.0   0.062
#frailty(ccode)                          16.4103 9.2   0.064

#Iterations: 7 outer, 48 Newton-Raphson
#Variance of random effect= 0.488   I-likelihood = -103.5 
#Degrees of freedom for terms= 0.7 0.8 0.7 0.8 0.8 0.8 0.7 0.7 9.2 
#Likelihood ratio test=105  on 15.1 df, p=1.78e-15
#n=597 (116 observations deleted due to missingness)

#T-G Test for proportionality 
cox.zph(cox.auto, transform='log')

#                rho chisq     p
#seq_ln       0.0726 0.298 0.585
##polity2    -0.1677 1.611 0.204
#public_cor   0.0491 0.203 0.652
#physical_vd -0.1659 1.418 0.234
#express_vd   0.1562 1.366 0.242
#intensity2   0.0639 0.204 0.652
#info         0.0858 0.540 0.462
#pop_ln      -0.1544 1.416 0.234
#GLOBAL           NA 6.600 0.580

#Graph Scoenfeld Residuals
ggcoxdiagnostics(cox.auto, type = "schoenfeld",
                 linear.predictions = FALSE,ggtheme = theme_bw())

#Result: model/individual variables do not violate proportionality assumption
#Diagnostics: Overall, model should be good because Global (p >.05). 

#ProportionalHazard
coef <- summary(cox.auto)$coefficients[1, 1]
coef
exp(coef)
#0.6012533
1-0.6012533
#0.3987467; ~40%

#########################
#####Anocracy############
#########################

#Estimate model
cox.ano <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                 + physical_vd + express_vd + intensity2 
                 + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.ano, 
                 control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.ano

#                   coef se(coef)      se2    Chisq DF       p
#seq_ln         -0.27733  0.10430  0.09772  7.06958  1 0.00784
#polity2         0.03201  0.03491  0.03155  0.84077  1 0.35918
#public_cor     -0.15536  0.64708  0.55041  0.05765  1 0.81026
#physical_vd    -3.06127  0.76923  0.66842 15.83757  1 6.9e-05
#express_vd      2.63870  0.77482  0.67584 11.59782  1 0.00066
#intensity2      0.80672  0.12868  0.11917 39.30533  1 3.6e-10
#info            0.01826  0.00788  0.00619  5.36749  1 0.02052
#pop_ln          0.30013  0.09914  0.08052  9.16537  1 0.00247
#frailty(ccode)                            43.10915 23 0.00663

#Iterations: 7 outer, 35 Newton-Raphson
#Variance of random effect= 0.348   I-likelihood = -425.2 
#Degrees of freedom for terms=  0.9  0.8  0.7  0.8  0.8  0.9  0.6  0.7 23.0 
#Likelihood ratio test=253  on 29 df, p=0
#n=1067 (107 observations deleted due to missingness)
#plot(survfit(res.cox),xlim=c(1992, 2014))

#T-G Test for proportionality 
cox.zph(cox.ano, transform='log')

#                 rho   chisq      p
#seq_ln        0.06343 0.69956 0.4029
#polity2      -0.02158 0.08774 0.7671
#public_cor    0.01390 0.04343 0.8349
##physical_vd -0.06030 0.83672 0.3603
#express_vd    0.00476 0.00478 0.9449
#intensity2   -0.16433 4.68451 0.0304
#info          0.02080 0.10811 0.7423
#pop_ln        0.02787 0.19644 0.6576
#GLOBAL             NA 6.77781 0.5608

#Graph Scoenfeld Residuals
ggcoxdiagnostics(cox.ano, type = "schoenfeld",
                 linear.predictions = FALSE,ggtheme = theme_bw())

#Result: intensity2 does not violate proportionality assumption
#Diagnostics: Overall, model should be good because Global (p >.05). 

#Proportional Hazard
coef <- summary(cox.ano)$coefficients[1, 1]
coef
exp(coef)
#0.7578026
1-0.7578026
#0.2421974; ~24%

##########################
#####Democracy############
##########################

#Estimate model
cox.demo <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.demo, 
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.demo

#                   coef se(coef)      se2    Chisq   DF       p
#seq_ln          0.13357  0.10871  0.09308  1.50981  1.0  0.2192
#polity2        -0.26551  0.09918  0.08649  7.16660  1.0  0.0074
#public_cor      1.22361  0.61004  0.50007  4.02321  1.0  0.0449
#physical_vd    -4.91608  0.90391  0.75152 29.57911  1.0 5.4e-08
#express_vd      2.49755  1.30963  1.11491  3.63688  1.0  0.0565
#intensity2      0.45261  0.15114  0.12094  8.96820  1.0  0.0027
#info            0.02479  0.00808  0.00618  9.41559  1.0  0.0022
#pop_ln          0.55187  0.08277  0.06538 44.46091  1.0 2.6e-11
#frailty(ccode)                            43.82388 23.5  0.0066

#Iterations: 7 outer, 66 Newton-Raphson
#Variance of random effect= 0.331   I-likelihood = -626.6 
#Degrees of freedom for terms=  0.7  0.8  0.7  0.7  0.7  0.6  0.6  0.6 23.5 
#Likelihood ratio test=565  on 28.9 df, p=0
#n=1922 (247 observations deleted due to missingness)
#plot(survfit(res.cox),xlim=c(1992, 2014))

#T-G Test for proportionality 
cox.zph(cox.demo, transform='log')

#                  rho    chisq       p
#seq_ln      -0.076110 2.22e+00 0.13646
#polity2     -0.138337 6.93e+00 0.00847
#public_cor  -0.095218 3.89e+00 0.04863
#physical_vd -0.011067 5.35e-02 0.81709
#express_vd   0.040459 6.40e-01 0.42359
#intensity2  -0.042492 8.74e-01 0.34979
#info        -0.056509 1.57e+00 0.20977
#pop_ln      -0.000682 2.08e-04 0.98849
#GLOBAL             NA 1.71e+01 0.02878

#Graph Scoenfeld Residuals
ggcoxdiagnostics(cox.demo, type = "schoenfeld",
                 linear.predictions = FALSE,ggtheme = theme_bw())

#Result: Model violates proportionality assumption (GLOBAL p<.05). Also, polity2 and public_cor.
#Diagnostics: Address polity2 and public_cor using spline regressions (Keele 2010).

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
######### SPLINE THERAPY ########
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

#R crashes, so run sperate spline regressions on polity, then public_cor to inspect results
# then run w/ the interactions for manuscript model

#Run nonlinearity tests 
cox.demo.spline.polity <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + pspline(polity2, df=3)
                                + public_cor
                                + physical_vd + express_vd 
                                + intensity2 
                                + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.demo, 
                                control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

summary(cox.demo.spline.polity)

#Public Cor
cox.demo.spline.pc<- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 
                           + pspline(public_cor, df=2)
                           + physical_vd + express_vd 
                           + intensity2 
                           + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.demo, 
                           control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

summary(cox.demo.spline.pc)

###Run new model w/ time interactions

#create time interactions
t<-sa.one.demo$year
corpt<-sa.one.demo$public_cor
polity<-sa.one.demo$polity2
sa.one.demo$yearxcorpt <- log(t)*corpt
sa.one.demo$yearxpolity <- log(t)*polity

cox.demo.spline2 <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                          + physical_vd + express_vd  + intensity2
                          + info + pop_ln + yearxcorpt + yearxpolity + frailty(ccode), data =  sa.one.demo,
                          control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.demo.spline2

#                   coef  se(coef)       se2     Chisq   DF       p
#seq_ln          1.43e-01  1.10e-01  9.31e-02  1.69e+00  1.0  0.1935
#polity2         4.78e+02  1.78e+02  1.73e+02  7.24e+00  1.0  0.0071
#public_cor      1.24e+03  8.86e+02  8.53e+02  1.96e+00  1.0  0.1612
#physical_vd    -4.90e+00  9.07e-01  7.52e-01  2.92e+01  1.0 6.7e-08
#express_vd      2.57e+00  1.32e+00  1.12e+00  3.76e+00  1.0  0.0526
#intensity2      4.72e-01  1.55e-01  1.24e-01  9.23e+00  1.0  0.0024
#info            2.48e-02  8.47e-03  6.22e-03  8.54e+00  1.0  0.0035
#pop_ln          5.62e-01  8.42e-02  6.64e-02  4.45e+01  1.0 2.6e-11
##yearxcorpt     -1.63e+02  1.16e+02  1.12e+02  1.96e+00  1.0  0.1615
#yearxpolity    -6.29e+01  2.34e+01  2.28e+01  7.25e+00  1.0  0.0071
#frailty(ccode)                                4.58e+01 24.2  0.0051

#Iterations: 7 outer, 56 Newton-Raphson
#Variance of random effect= 0.355   I-likelihood = -622.8 
#Degrees of freedom for terms=  0.7  1.0  0.9  0.7  0.7  0.6  0.5  0.6  0.9  1.0 24.2 
#Likelihood ratio test=575  on 31.9 df, p=0
#n=1922 (247 observations deleted due to missingness)

#T-G Test for proportionality 
cox.zph(cox.demo.spline2, transform='log')

#                rho  chisq      p
#seq_ln      -0.0977  3.789 0.0516
#polity2     -0.1240  2.874 0.0900
#public_cor  -0.0554  0.669 0.4132
#physical_vd  0.0224  0.217 0.6411
#express_vd   0.0200  0.159 0.6899
#intensity2  -0.0235  0.276 0.5993
#info        -0.0714  2.729 0.0985
#pop_ln       0.0306  0.431 0.5113
#yearxcorpt   0.0554  0.668 0.4136
#yearxpolity  0.1240  2.875 0.0900
#GLOBAL           NA 11.387 0.3281

#Graph Scoenfeld Residuals
ggcoxdiagnostics(cox.demo.spline2, type = "schoenfeld",
                 linear.predictions = FALSE,ggtheme = theme_bw())

#Result: After using splines, model/variables do not violate proportionality assumption.
#Diagnostics: Overall, model should be good because Global (p >.05). 

#Create data for Table 10, Appendix#

#########
###ONE###
#########
rm(list=ls())

#Load data
sa.one.full<-read.csv("survival_one.csv") 

#Stratify sample

#Autocracy
sa.one.auto <- subset(sa.one.full, durable2 == 0)
#Anocracy
sa.one.ano <- subset(sa.one.full, durable2 == 1)
#Democracy
sa.one.demo <- subset(sa.one.full, durable2 == 2)

#Make duration data. Killed = events/failure
sa.one.full$killed <- with(sa.one.full, Surv(year, one == 1))

#How many events?

#Global
one<-sa.one.full$one
CrossTable(one) 
#439

#Auto
one<-sa.one.auto$one
CrossTable(one) 
#48

#Ano
one<-sa.one.ano$one
CrossTable(one) 
#154

#Demo
one<-sa.one.demo$one
CrossTable(one) 
#219

#########
###TWO###
#########

#Load data
sa.two.full<-read.csv("survival_two.csv") 

#Stratify sample

#Autocracy
sa.two.auto <- subset(sa.two.full, durable2 == 0)
#Anocracy
sa.two.ano <- subset(sa.two.full, durable2 == 1)
#Democracy
sa.two.demo <- subset(sa.two.full, durable2 == 2)

#Make duration data. Killed = events/failure
sa.two.full$killed <- with(sa.two.full, Surv(year, two == 1))

#How many events?

#Global
two<-sa.two.full$two
CrossTable(two) 
#203

#Auto
two<-sa.two.auto$two
CrossTable(two) 
#18

#Ano
two<-sa.two.ano$two
CrossTable(two) 
#68

#Demo
two<-sa.two.demo$two
CrossTable(two) 
#101

###########
###THREE###
###########

#Load data
sa.three.full<-read.csv("survival_three.csv") 

#Stratify sample

#Autocracy
sa.three.auto <- subset(sa.three.full, durable2 == 0)
#Anocracy
sa.three.ano <- subset(sa.three.full, durable2 == 1)
#Democracy
sa.three.demo <- subset(sa.three.full, durable2 == 2)

#Make duration data. Killed = events/failure
sa.three.full$killed <- with(sa.three.full, Surv(year, three == 1))

#How many events?

#Global
three<-sa.three.full$three
CrossTable(three) 
#119

#Auto
three<-sa.three.auto$three
CrossTable(three) 
#11

#Ano
three<-sa.three.ano$three
CrossTable(three) 
#47

#Demo
three<-sa.three.demo$three
CrossTable(three) 
#52

#Creat figures for Figure 6, appendix#

#Figure 6a-One journalist killed threshold

#Rerun models for figure
cox.auto.one <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.auto,      
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

cox.ano.one <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                 + physical_vd + express_vd + intensity2 
                 + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.one.ano, 
                 control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

#create time interactions
t<-sa.one.demo$year
corpt<-sa.one.demo$public_cor
polity<-sa.one.demo$polity2
sa.one.demo$yearxcorpt <- log(t)*corpt
sa.one.demo$yearxpolity <- log(t)*polity

cox.demo.one.spline2 <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                          + physical_vd + express_vd  + intensity2
                          + info + pop_ln + yearxcorpt + yearxpolity + frailty(ccode), data =  sa.one.demo,
                          control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))
####Plot three regime types
plot(survfit(cox.auto.one), col = c('black') , conf.int=F,lwd=4, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.ano.one), col = c('black'), conf.int=F, lty="dotted",xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.demo.one.spline2), col = c('black') , conf.int=F, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
legend(1986.5,0.35, c("Autocracies","Anocracies", "Democracies"),
       lty = c("solid","dotted","solid"),lwd=c(4,1,1), col=c("black","black", "black"), 
       cex=2.25, y.intersp = .25,x.intersp = .15,seg.len=.5, bty = "n", adj = c(0, .5))
dev.off(dev.list()["RStudioGD"])

#Figure 6b-Twojournalists killed threshold
cox.auto.two <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data = sa.two.auto,      
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))
cox.ano.two <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                 + physical_vd + express_vd + intensity2 
                 + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.two.ano, 
                 control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))
cox.demo.two <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =  sa.two.demo, 
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

####Plot three regime types: Two killings threshold
plot(survfit(cox.auto.two), col = c('black') , conf.int=F,lwd=4, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.ano.two), col = c('black'), conf.int=F, lty="dotted",xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.demo.two), col = c('black') , conf.int=F, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
legend(1986.5,0.35, c("Autocracies","Anocracies", "Democracies"),
       lty = c("solid","dotted","solid"),lwd=c(4,1,1), col=c("black","black", "black"), 
       cex=2.25, y.intersp = .25,x.intersp = .15,seg.len=.5, bty = "n", adj = c(0, .5))
dev.off(dev.list()["RStudioGD"])


#Figure 6c-Threejournalists killed threshold

#Autocracies gives warning, but like my Cajun uncle says, "It won't hurt none" 
cox.auto.three <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data =sa.three.auto,      
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))
cox.ano.three <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                 + physical_vd + express_vd + intensity2 
                 + info + pop_ln + frailty(ccode), na.action=na.exclude, data = sa.three.ano, 
                 control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))
cox.demo.three <- coxph(Surv(X_t0, X_t, X_d) ~ seq_ln + polity2 + public_cor
                  + physical_vd + express_vd + intensity2 
                  + info + pop_ln + frailty(ccode), na.action=na.exclude, data = sa.three.demo, 
                  control=coxph.control(eps=1e-09,iter.max=100,outer.max=100))

####Plot three regime types
plot(survfit(cox.auto.three), col = c('black') , conf.int=F,lwd=4, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.ano.three), col = c('black'), conf.int=F, lty="dotted",xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
par(new=TRUE)
plot(survfit(cox.demo.three), col = c('black') , conf.int=F, xlim=c(1990, 2015), ylab = "Survival Probability",
     xlab = "Years",cex.lab=1.75, mgp=c(2,.75,0))
legend(1986.5,0.35, c("Autocracies","Anocracies", "Democracies"),
       lty = c("solid","dotted","solid"),lwd=c(4,1,1), col=c("black","black", "black"), 
       cex=2.25, y.intersp = .25,x.intersp = .15,seg.len=.5, bty = "n", adj = c(0, .5))
dev.off(dev.list()["RStudioGD"])

################
##  END ~ jas ##
################